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obtain simple expressions for the dependence of the final particle number on the expansion velocity 
and final dielectric constant. 
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I. INTRODUCTION 

One of the most intriguing and fascinating features of 
quantum field theory resides in the non-trivial nature 
of its vacuum states. Quantum fluctuations present in 
the vacuum are responsible for non-classical effects that 
can be experimentally detected. The most well known of 
such phenomena is perhaps the Casimir effect [1]. In this 
setting, the vacuum energy between two static metallic 
plates is shifted due to quantum fluctuations, resulting 
in an attractive force between the plates. In general non- 
static situations, changes in a system's parameters may 
lead to a redefinition of the natural vacuum state. As 
a result, a system initially in the vacuum may at later 
stages display a non-zero particle content. The possi- 
bility of this type of particle production mechanism has 
been identified and studied in a large variety of systems. 
In a cosmological context, particle creation may occur 
in back- hole radiation - the Hawking effect [2], and as 
a consequence of an expanding gravitational background 
[3]. In the laboratory it is possible that dynamical gen- 
eralisations of the Casimir effect may lead to observable 
effects [4]. Alternatively, there have been efforts to re- 
produce the Hawking effect using specific matter systems 
[5,6]. This is the case of slow-light experiments [7], black- 
hole sonic analogues in Bose-Einstein condensates [8] and 
dielectric media [9]. 

The calculations involved in determining the physical 
outcome of particle creation processes, though trivial to 
state, are often hard or impossible to complete. Usu- 
ally one is required to find the solution of a set of space 
and/or time-dependent field equations, with initial condi- 
tions covering a complete basis of functions. Even when 
relying on simplifying approximations, the set of prob- 
lems for which solutions can be found is considerably 
limited. On the other hand, linear partial differential 
equations are relatively easy to solve numerically, and the 
above task is within present computational capabilities. 



In this paper we explore this possibility and introduce a 
fully numerical scheme for studying particle production 
in general settings. We will use examples with known 
solutions to help developing and testing the required nu- 
merical techniques. These will then be applied to new 
cases, extending original predictions and leading to new 
results. 

In Section II, we rely on a simple cosmological exam- 
ple to introduce notation and discuss concepts relevant 
for the rest of the paper. We use this system to describe 
the general numerical approach and we test it, by com- 
paring the results with known analytical predictions. In 
Section III, we apply these techniques to the case of par- 
ticle creation by a reflecting moving boundary in a one- 
dimensional cavity. For an oscillating motion we obtain 
particle production rates and spectra, for both resonant 
and non-resonant frequencies, in the short and long-time 
regimes. In Section IV, we study a system with time 
dependent optical properties, mimicking a bubble of di- 
electric material expanding into the vacuum. From the 
numerical data we determine the dependence of the final 
particle number on the expansion velocity and strength of 
the dielectric. Finally, in Section V we discuss our results 
and suggest further applications. Appendix A contains a 
description of the numerical algorithm used in Section III 
to solve the field equations for moving boundary condi- 
tions. 

II. A FIRST EXAMPLE: PARTICLE CREATION 
IN AN EXPANDING BACKGROUND 

In this section we will look at a simple example of 
particle creation in a cosmological system. In particular, 
we will study the evolution of a free complex field <fi in an 
expanding flat universe. Here and throughout the paper 
we will follow the conventions and notation of reference 
[3]. 
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For a general gravitational background, the equation 
of motion for a scalar field is given by 



(-g)- 1/2 d^(-g) 1/2 g^d u( f>} + m 2 <p = 0, 



(1) 



where g M „ is the space-time metric and g — | det g„ v \ . 

The invariant scalar product of any two solutions of 
Eq. (1) is defined as 



hdpfo-d^KfcW^dE' 1 , (2) 



where the integral is taken over any space-like surface 
of the space-time is the time-like unit-vector or- 

thogonal to the surface). Using Eq. (2) we can construct 
orthonormal basis of solutions of the equations of motion, 
obeying: 

(u i ,u j )=6 ij , {u*,u*) = -6ij, (ui,u*)=0 (3) 

Expanding the field operator in terms of a given basis we 
obtain the familiar creation and annihilation operators 



0( x ) = E ^ u ^ x ) + a l u i( x )] 



(4) 



Particle states are defined in respect to a\, aj. In general 
we will be looking at systems which are static and flat 
for both early and late times - the so called in-out type 
problems. For these cases we can single-out in a natu- 
ral way, two different basis of solutions, {uf 1 } and {w° ut } 
which tend respectively as t — > U n and t — > t out , to the 
usual flat-space positive frequency modes. The field can 
be quantised in terms of any of these two basis, leading 
to different definitions of creation and annihilation oper- 
ators and hence to different definitions of particle states. 
In most cases of interest, the particle content of a given 
state will change depending on the set of modes we chose 
as basis. In particular, we will consider the case where 
the system is in the vacuum state corresponding to the 
in modes. In respect to the out particle states - the natu- 
ral choice for late times - this state may have a non-zero 
particle content. The particle spectrum is given by [3] 



(5) 



where iV, is the number of particles in the mode i. The 
(3ij 's are a subset of the so-called Bogolugov coefficients 
which relate the two basis and can be evaluated using [3] : 



a. 



= « Ut ,<), 



p ij = -(ur t ,[<Y 



(6) 



In this paper, we will look into how to tackle the cal- 
culation of the particle spectrum for general in-out prob- 
lems, using a fully numerical approach. We will start 
by evaluating numerically the in modes {u- n }. This is 



done by solving the equation of motion Eq. (1), taking 
as initial condition the asymptotic, flat-space mode so- 
lutions corresponding to the Minkowsky vacuum. These 
are propagated up to t out , when the new stationary flat 
configuration is reached. The resulting in modes can then 
be contracted with the out modes, by evaluating numer- 
ically the integral Eq. (2) at t — t out - In this way we 
obtain the Bogolugov coefficients Eq. (6). Other physical 
quantities, such as the particle spectrum can be evaluated 
in terms of these. 

To illustrate the method, we look at a simple example 
of a 1 + 1 space-time with a time dependent scale factor. 
In conformal coordinates, the metric is given by 



d.s 2 =a 2 (r 1 )(dr 1 2 -dx 2 ), 
and the field equations of motion are 

d 2 (j> d 2 ^ 



drf dx 2 



— ma {rj)4> 



(7) 



(8) 



For this system, the time independent inner product of 
two solutions is given by 



(0i 



/+oo 
-OO 



dx 



(9) 



We take a{rf) — > a m as rj — > — oo, and a(rf) — > a ou t as 
r\ — > +oo. In these limits, the corresponding mode basis 
have the usual flat space-time form: 



1 



1 



ikx—iuj\ n r) 



+ 171 (I 



1 n 1 



(10) 



ikx-iuioutV 



2y/lTLd ou t 



W ou t 



+ m 2 a 2 out (11) 



For the scale- factor a(i]) we choose a particular form for 
which there is a known analytical solution [3] : 



a 2 ( V ) 



+ afn) + o ("out - «in) tanh(p?7) (12) 



Taking Eq. (10) as initial condition, we evolve the equa- 
tion of motion Eq. (8) up to a time when a ~ a out . 
The numerical solutions are then contracted with the out 
modes Eq. (11), using Eq. (9) to obtain the Bogolugov 
coefficients. For our choice of scale-factor these can be 
evaluated analytically. The derivation can be found in 
[3], the result being given by 



[awl 2 = Skk> 



>kk> 



sinh 2 [7r(o; out +w in )/(2p)] 
sinh(7rcj in //9) sinh(7rw out /p) 

sinh 2 [7r(u out - m in )/(2p)} 
sinh(7rw in /p) sinh(7rw ou t//3) 



(13) 
(14) 



In Fig. 1 we compare the numerical Bogolugov coef- 
ficients with the expected analytical values, for a par- 
ticular choice of parameters m = 1, p = 10., di n = 1. 
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and a out = 10. We evolved 40 modes for a period of time 
iout — tin = 0.6, enough for the scale factor to change from 
its initial to final value with good accuracy. The equation 
of motion Eq. (8) was discretized using a standard Leap- 
Frog algorithm, in a box of physical size L = 100 with 
periodic boundary conditions. The lattice spacing and 
time-step were set to dx = 0.25 and dt = 0.0001. The 
numerical data matches very well the analytical predic- 
tion, confirming the accuracy of the numerical approach 
for the chosen set of parameters. 
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FIG. 1. Analytical (line) and numerical results (crosses) 
for Bogolugov coefficients, for a model with m = 1, p = 10., 
ai n = 1. and a out = 10. 

Note that only diagonal elements of the Bogolugov ma- 
trix are non-zero, a consequence of the invariance of the 
system under spatial translations. This implies that the 
particle spectrum coincides with the diagonal of the Bo- 
golugov matrix, i.e. Nk = |/3fefc| 2 - Since Eq. (8) is sep- 
arable in space and time, we would have the choice of 
evolving independently each mode's amplitude in Fourier 
space. Numerically this approach would have been con- 
siderably simpler. Our method is crucial though, in sit- 
uations where this option is not available. In the next 
sections we will apply it to a few of such cases, systems 
that have both explicit space and time-dependency. 



III. PARTICLE PRODUCTION BY A MOVING 
MIRROR IN A CAVITY 

In this section we will look at particle production by 
a moving mirror in a closed cavity. This system has at- 
tracted the attention of several authors, displaying many 
complex and non-trivial features despite its apparent sim- 
plicity. More importantly, it is the most natural setting 
for an experimental observation of particle creation from 
a vacuum state (see [4] and references therein). 

We will consider a scalar field confined to a one- 
dimensional box with one of the boundaries fixed at x = 



and the other moving according to a given trajectory 
X(t) (with \X{t)\ < 1). For the sake of simplicity, and 
in order to take advantage of conformal invariance we 
assume the field to be massless, its equation of motion 
being given by 



d 2 <t> 
dt 2 



d 2 <f> 

dx 2 



(15) 



The boundary conditions for the field are obtained by 
constraining it to vanish on both walls at all times, i.e. 
0(0, t) — (j)(X(t), t) — 0. By doing so we assume the walls 
to be perfect reflectors. 

In the case of a stationary mirror, time translational 
invariance makes it possible to find a set of mode solu- 
tions with positive definite frequency. These are given 

by 



u„(x,t) 



1 



n-K 



; sin(w„x), 



rnr 



ui n = — n = 1, ... ,oo 



(16) 



where L is the length of the static box. These modes 
and their complex conjugates form a complete set of solu- 
tions of Eq. (15), orthonormal in respect to the simplectic 
product 



(01,^2) 



• L 



[01 <9 t 02 - d, 



dx 



(17) 



The stationary solution allows us to define a general 
class of in- out problems. We will consider mirror trajec- 
tories such that X{t) is a constant X for t < 0, evolves 
arbitrarily for a period of time, going back to rest at Xf 
for t > tf. The time and space dependence introduced in 
the system by the motion of the mirror, leads to particle 
production without the need to consider additional exter- 
nal fields. This phenomenon is known as the dynamical 
Casimir effect. 

In order to determine the final particle content of the 
box for a given X(t), one must find a set of solutions for 
the equation of motion, using Eq. (16) as initial con- 
dition. This problem can be formally solved by tak- 
ing advantage of the conformal invariance of Eq. (15). 
G. Moore showed [10] (see also [11]) that by changing 
coordinates, the general problem can be mapped into 
the stationary case. In this way a complete set of so- 
lutions can be found, with the general form u n (x,t) cx 
cxp[— inirR(t — x)\ — exp[—innR(t + x)]. R{u) is a func- 
tion to be determined in terms of the particular motion 
of the mirror, according to [10] 



R[t + X(t)] = R[t-X(t)} + 2 



(18) 



This equation is of course non-trivial to solve in gener- 
ality. One approach is to find R(u) in terms of X(t) via 
a perturbative expansion in the velocity of the mirror 
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[10]. This technique is valid only in the adiabatic limit 
X « 1 and gives no meaningful results in the relativis- 
tic regime. Only for a few specific cases of the mirrors's 
trajectory have analytical solutions been found [12-15]. 

Using the numerical procedure outlined in Section II 
the problem can be solved for any given arbitrary mirror 
trajectory. The known analytical solutions will be useful 
in providing tests for the algorithms and in gauging the 
numerical requirements for good accuracy. The general 
procedure will be as before: the m-modes Eq. (16) are 
evolved numerically up to the final time tf, when the re- 
sult is contracted with the basis of stationary modes for 
the final length Xf. The main difficulty lies in obtaining 
a numerical solution for the wave equation with moving 
boundaries. Since the size of the box changes with time 
we are forced to vary the number of points in the simula- 
tion lattice. These discontinuous jumps in the lattice size 
introduce errors in the numerical solution. We overcame 
the problem by developing an improved Leap-Frog algo- 
rithm that smoothes the values of the field derivatives on 
the boundary - a detailed discussion of the method can 
be found in Appendix A. 

Finally we should stress that although the focus here 
will be on one dimensional massless field theories, the 
same methods can easily be extended to higher dimen- 
sions and massive fields [16]. In those cases very few re- 
sults are known [17] since one cannot rely on conformal 
invariance to simplify the problem. For those systems, a 
numerical approach may be the only way to tackle accu- 
rately the long time regime. 



A. Uniformly Moving Mirror 

Here we will look briefly at a simple type of trajectory 
for which the in-out problem can be solved analytically 
[10-12], namely the case when the mirror moves with 
constant velocity. The mirror evolves from rest at t = 0. 
moving according to X(t) = Xq + vt up to final time tf, 
when it becomes stationary again. The function R(u) 
obeying Eq. (18) for this trajectory was fully determined 
by Castagnino and Ferraro [12]. These authors showed 
that R is a piece-wise linear function with derivative dis- 
continuities. This is a consequence of the discontinuous 
change in the velocity of the mirror for initial and final 
times. The final spectrum for created particles was es- 
timated to be of the form N(k) oc v 2 /k, in the limit of 
large k and velocities much lower than the speed of light. 

In Fig. 2 we confirm this estimate. The results were 
obtained by simulating a box with initial physical length 
Xq = 50. expanding with constant speed v for a total 
time tf — 50. We evaluated the final particle spectrum 
for 20 modes, for values of v ranging from 0.1 to 0.9. As 
seen in the figure, the spectrum for high k is well fitted 
by a power law of the form N(k) — Bk a , with B and a 
depending on v. As expected, a is close to —1, varying 



between —1.1 for v = 0.1 and —1.3 for v = 0.9. B also de- 
pends quadratically on the expansion velocity for small v. 
When relativistic speeds are reached (about v <~ .7) this 
simple relation stops holding, with B increasing faster 
and leading to a stronger rate of particle production. 

Note that the 1/k spectrum leads to a logarithmically 
divergent total particle number. This divergence is a con- 
sequence of the discontinuous changes in the velocity of 
the mirror [10]. As a rule of thumb we may expect parti- 
cle production to be related to the mirror's acceleration. 
As this diverges for the initial and final times, so docs 
the total particle number. In realistic cases however, the 
divergence should disappear, as the acceleration remains 
bounded for all times. 
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FIG. 2. Left plot: particle spectrum for mirror velocities 
varying from 0.1 to 0.9 (bottom to top). The straight lines are 
fits to a power law N = Bk a , using the 15 higher momentum 
data points. In the right plot we have the dependence of B 
on the velocity. For low-i> this is well described by B = v 2A , 
the numerical power being obtained by fitting the first 6 data 
points. 

Finally let us note that for a box that expands contin- 
uously with constant v, with no initial or final stationary 
points, an exact solution for Eq. (18) can easily obtained 
[10]: 



R{u) = 



arctan v 



log [u + 



Xc 



(19) 



with X the position of the mirror at t = 0. The corre- 
sponding normalised mode basis is given by 



2-^/wr 



-iQ n \og(Xo/v+t— x) 



ifl n log(X /v+t+x) 



nir 



arctanh(u) 



n = 1, 



(20) 



This set of exact analytical solutions proved useful in 
testing the accuracy of the numerical algorithm for the 
moving boundary problem, as detailed in Appendix A. 
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An interesting feature of the solution is that for t r = 
2Xo/(l — v) the evolved modes coincide with the origi- 
nal ones, re-scaled to the new box length. This feature 
is also present for the in-out motion studied above. It 
implies that at t r the in modes are identical to the out 
modes and the Bogolugov coefficients should vanish. As 
a consequence we should observe no particle production 
for this particular stopping time. We checked that this 
'resonance' does indeed take place by performing the sim- 
ulation described above with final time set to tf = t r . 
For several choices of expansion velocity (since t r grows 
with the inverse of 1 — v leading to very long simulation 
times, we restricted ourselves to the \ow-v regime) the 
total number of particles produced was indeed zero. 

B. Oscillating Mirror 

Moving a step up in complexity from the uniformly 
moving mirror problem, we will consider in this section 
a cavity whose boundary executes small, periodic oscilla- 
tions. This type of dynamical Casimir effect setting has 
been widely studied, and there is hope that an equivalent 
of a vibrating mirror could be experimentally realized. In 
the particular situation where the mirror frequency res- 
onates with one of the cavity modes, particle production 
effects may be magnified, leading eventually to observ- 
able results. A thorough review of the topic with exten- 
sive references can be found in [4]. 

We will consider a sinusoidal mirror trajectory given 

by 

X(t) =X + ^[1- cos(wt)] (21) 

In the simulations described below, the amplitude of the 
motion A will be taken to be considerably smaller than 
Xq, the initial size of the box. Typically we will set 
Xq = 50 and A — 2. The frequency of the motion uj, 
will vary for different runs, care being taken that the ve- 
locity of the mirror never exceeds the speed of light. In 
general to will be set to be one of the natural frequencies 
of the cavity, ui n — mr/XQ. This choice of frequencies 
is usually favoured in the literature, since it is expected 
that resonance between the cavity modes and the motion 
of the mirror will lead to higher particle production for 
long times. The resonant case is also considerably easier 
to tackle analytically (see [4] where the problem is solved 
exactly for n = 2 and [18] for a calculation of R(u) for 
a family of periodic-like trajectories). In fact, as far as 
we are aware, off-resonant solutions are only discussed 
in [19], where the energy production is evaluated for the 
general case. Our method of course, does not distinguish 
between resonant and off-resonant frequencies - in gen- 
eral we will present results for resonant frequencies and 
discuss how they change for the general case. 



We start by looking at results for very short times, 
when the mirror executes a small number of oscillations. 
In Fig. 3 we show the total number of particles produced 
as a function of time, for several choices of the oscillation 
frequency. At every time-step, the numerical field profiles 
were contracted with the static basis for a stationary box 
of length X(t). As discussed in the previous section this 
will lead to a convergent particle spectrum, but only for 
times when the mirror velocity is zero. For the sake of 
clarity, and with this caveat in mind, we still plot the 
total particle number for all time steps. For the number 
of modes simulated, the results for non-stationary final 
states do not deviate much from the finite case. Since 
the divergence is logarithmical, we presume it would only 
manifest itself if much higher frequencies were included. 
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FIG. 3. Total number of particles versus time for different 
frequencies of the oscillating mirror. The straight lines are 
linear fits to iVtot (t) . The oscillation frequencies are of the 
form ui = nir/Xo with n — 4,6, 8, 10, 12, 14. 60 modes were 
used to calculate the total particle number. 

It is clear from Fig. 3 that the evolution is effectively 
linear for the period shown. The production rate changes 
only for times of the order of twice the box length - up 
to then each oscillation produces roughly the same in- 
crease in particle number. This suggests that for small 
times, the field disturbances radiate away from the mov- 
ing boundary without significant interference. For later 
times, after reflecting back on the stationary wall, they 
will interact with the vibrating mirror, constructively 
or otherwise depending on their frequency, leading to a 
change in the production rate. 

For each value of the oscillating frequency u>, we fitted 
the results to a linear law of the form AT tot (i) = jt. We 
checked that allowing for a constant term in the linear 
expression, AT to t (t) = b + -ft did not change the estimate 
for the production rate 7 (b was always found to be very 
small, of order less than 10~ 3 ). Clearly, higher oscillating 
frequencies lead to higher particle production rates. In 
Fig. 4 the dependence of 7 on the mirror's frequency is 
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made more explicit. The results are very well fitted by 
a power law of the form N tot (t) — Du k t with k = 3.3, 
a result in agreement with the analytical estimates in 
the literature. It is generally predicted that for a single 
oscillation the number of particles produced varies as a 
power law of the oscillation frequency. The particular 
value of the power depends on the specific type of oscil- 
lating trajectory. In [12] the production power for half 
an oscillation was shown to vary between 2. and 4., and 
in [13] a u> 4 dependence was found in the limit of small 
velocities. For a single oscillation, our result corresponds 
to a power of 2.3, which is in line with these values. 
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FIG. 4. Fit of the particle production rate to a power law, 
log- log and linear scales. 
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FIG. 5. Total particle production versus time for the n=2 
resonance frequency. The fit to a power law - dashed line - is 
for times varying from t = 100 till t = 700. 40 modes were 
used in the simulation. 

As for the amplitude coefficient D, it was suggested in 
[13] that it should vary quadratically with the amplitude 
of the vibration of the mirror. We performed a quick 
test of this estimate by repeating the same set of runs 



with A, the mirror oscillating amplitude set to 1., half its 
previous value. Fitting both sets of data simultaneously 
to the same power k and allowing D to vary, we obtained 
k = 3.2, D = 0.0342 and D = 0.00780 for A = 2. and 
A = 1. respectively. The rate of the two amplitudes 
is 0.23, consistent with a reduction by a factor of 4 as 
expected. 

We now look at the longer time behaviour of the sys- 
tem. As mentioned above, the evolution stops being de- 
scribed by a linear law at about t ~ 2A , when a new 
regime with a different production rate sets in. In Fig. 5 
we can see the total particle number for longer times for 
the resonant frequency u> = 2ir / 'X$. 

The evolution for the period of time shown is well fitted 
by a power law JV tot = 8.5 x 10~ 7 i", with a = 1.85. For 
the dynamical range in question a power law growth with 
a = 2. was predicted in ref. [4]. This is in reasonable 
agreement with our numerical estimate. 

For higher frequencies, as far as we are aware, there 
are no explicit predictions for the particle number time 
dependency (see [18,19] for a calculations of the energy 
produced). In Fig. 6 we show N tot as a function of time 
for two higher resonant frequencies u) = 6tt/X and u = 
7n/Xo. We also include iV t ot f° r the off-resonant case 
u> = 6.5 tt/Xq. 
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FIG. 6. Total number of particles produced for resonant 
frequencies n=6 and n=7 and off-resonant frequency n=6.5. 
40 modes used in the simulation. 

The overall evolution does not fit any simple expres- 
sion in any of the three cases. We note that the initial 
period of linear evolution is always present and we have 
checked that the production rates follow the power law 
shown in Fig. 4, even in the off-resonant case. For later 
times though, the contrast between off-resonant and reso- 
nant trajectories is very clear, with limited quasi-periodic 
production in the former case and continuing increase in 
the later. It is interesting not to find any qualitative dif- 
ference between the long time behaviour of the even and 
odd resonant frequencies systems, n=6 and n=7. In fact 
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it has been often suggested that the main phenomena 
governing the evolution in the long time limit is akin to 
parametric resonance [4]. If that were strictly the case 
though, resonant amplification would only take place for 
oscillating frequencies equal to twice the cavity's natural 
frequencies [20]. That is, using our notation, we should 
see long time suppression of the production rate for odd 
values of n, which clearly is not the case. This seems to 
confirm the results in [19] where it was found that the 
energy produced grows exponentially for long times, for 
any motion with uj equal to a multiple of the fundamental 
frequency tt/Xq. 




FIG. 7. Particle spectrum for w = 6ir/X , t = 720. 

We now look at the distribution of the particles pro- 
duced in terms of their frequency. In Fig. 7 we show the 
late time particle spectrum for the resonant n = 6 mo- 
tion. We observe an alternating succession of excited and 
suppressed frequencies. In particular we see that there 
is no production at all for frequencies multiple of the 
oscillating frequency This cancellation is exact and it 
coincides strictly with points of the form Nw, with N in- 
teger. The peaks of the spectrum follow a similar but less 
precise pattern. The most excited frequency is uj/2 and 
the following peaks take place at about points Mui + ui/2, 
though there seems to be a tendency for the peak to move 
up for higher frequencies. The fact that the maximum of 
the production rate corresponds to a frequency of about 
half the excitation, with higher even multiples suppressed 
and odd ones enhanced, is once again reminiscent of para- 
metric resonance. Nevertheless, since the same pattern 
is observed for oscillating frequencies odd multiples of 
ir/Xf), we must conclude that a more complex process is 
operating in the system. 

Finally we contrast the resonant spectrum with the 
n = 6.5 off-rcsonant case as shown in Fig. 8. Clearly, the 
overall amplitude is greatly reduced (by about one order 
of magnitude) and the pattern of alternating peaks is no 
longer present. The single peak is once again centred 



around ui/2 but its height oscillates widely in time, lead- 
ing to the periodic fluctuations in total particle number 
shown in Fig. 6. 
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FIG. 8. Particle spectrum for the off-resonant motion 
uj = 6.57t/X , for t = 720. 



IV. NON-UNIFORM TIME- VARYING 
DIELECTRIC 

In the examples discussed above, the varying geometry 
of the systems considered was the primary cause behind 
particle production from the initial vacuum state. In this 
Section we will look at a system with fixed dimensions, 
but whose bulk physical properties change both in space 
and in time. In particular we will consider wave propa- 
gation in a background with a space and time dependent 
dielectric "constant" . The system will be modelled by a 
straightforward generalisation of the wave equation [21] 

d ( . ,0<A d 2 <j) 

where e(x,t) is the varying dielectric term. Other gener- 
alisations are possible, their form depending ultimately 
on the properties of the particular system one chooses to 
model. In [22] , the change in e is caused by a moving half- 
infinite dielectric material, and compliance with Lorentz 
invariance leads to a wave equation with a cross deriva- 
tive term. In so-called slow light experiments [7], based 
on electromagnetically-induced transparency, the optical 
properties of the medium are controlled by an external 
source. These systems can be described by an effective 
Lagrangian that includes a mass-like term for the field, 
leading to yet another choice of field theory. Here, for 
the sake of simplicity we follow [21] and use Eq. (22), 
though our approach can be easily extended to the cases 
mentioned above. 

For two solutions of Eq. (22) the time invariant scalar 
product is given by 
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t(x,t)[<h.a t <&-d t <h.<fo\dx (23) 



the fields being denned in a finite box of size L. We will 
look at cases where e(x, t) is spatially uniform both for 
the initial and final times, t\ and t 2 . The corresponding 
normalised mode basis are given by: 



u n (x,t) = \l —e lu)nt sin(fc„x), 

TVK 



n = 1, . . . , oo 



(24) 



where c, = l/^/el with i = 1,2 denotes the initial and 
final speed of light respectively. 

In the specific scenario we will be simulating, the 
original medium with dielectric constant e\ is replaced 
throughout the box by a medium characterised by e 2 . 
Medium 1 expands into medium 2 from left to right with 
a given velocity v, as a one-dimensional analogue of an 
expanding bubble. At intermediate time steps, the two 
regions with different dielectric constants are separated 
by a "wall" that interpolates between ei and e 2 . We 
model the wall profile as a sine 

{ei, if x > d/2 

e+ - e_ sin(7rx/d), if -d/2 < x < d/2 (25) 
62, if x < -d/2 

with e + = (ei + 62)/2, e_ = (ei — 62)/2, d being the 
wall thickness. Using the wall profile x{ x ) we can easily 
define a time-dependent dielectric function modelling the 
setting described above, as e(x,t) = x( x — y t) with t 
varying between t\ = —d/(2v) and t 2 — d/(2v) + L/v. 
We will be looking at cases where the final speed of light is 
less than the original one, corresponding to a vacuum box 
being filled with a denser medium. The initial dielectric 
constant will be set to e\ = 1 with e 2 > e\. In general we 
will constrain the wall expansion velocity to be less than 
the speed of light in both the initial and final medium, 
that is v < 1 / \f^2 < 1 / v 7 ^! 

In Fig. 9 we illustrate the dependence of the number 
of particles produced on the expansion velocity. The wall 
thickness was set to d = 5., considerably smaller than the 
physical box size, L — 50. The initial dielectric constant 
was fixed to ei = 1. whereas its final value e 2 was varied 
between 1.2 and 4.0. For each particular choice of e 2 we 
simulated a series of expansion velocities ranging from 
v ~ 0.3 up to the limit v ~ l/^/ez. 

We found that for low velocities the production rate 
was considerably suppressed. As the expansion velocity 
increases, the total number of particles produced goes up 
steeply - on average by more than two orders of magni- 
tude - as shown in Fig. 9. For mid and high values of the 
velocity, the dependence of the final particle number on 
v is well described by an exponential law, N to t — Ae 1v . 
This expression is valid up to values of v considerably 
near the speed of light in the second medium. 
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FIG. 9. Left plot: particle number vs expansion velocity 
for values of £2 varying between 1.2 and 4.0 (bottom to top 
curves). The data was fitted for medium and high-w to an 
exponential Ae lv . The dependence of 7 on the final dielectric 
constant £2 is shown in the right plot. The results are well 



fitted by a linear function 7 
and K\ =5.1 




^ + with K = 12.0 
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v d 
FIG. 10. Left plot: particle number dependence on the ex- 
pansion velocity, for three different values of the wall thickness 
d = 2, 5, 12 (top to bottom). For the two lower values of d we 
fitted the data to AT tot = Ae ,v . For the higher value of the 
thickness, the results no longer follow a simple exponential 
law. In the right plot we show 7 as a function of the wall 
thickness for d varying between 1 and 7. 

When the final dielectric constant is allowed to vary we 
observe as expected, that stronger changes in the optical 
properties of the system lead to higher production rates. 
As a consequence we should observe an increase in the 
coefficient 7 for higher values of e 2 . This is fact the case, 
and it turns out that the dependence of 7 on the final 
dielectric constant is very simple. In the regime studied 
7 varies linearly with e 2 as can be seen in the fit shown in 
the right plot of Fig. 9. Using dimensionality arguments 
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we can then write 



iVtot oc exp 



Kn 



£1/ Ci 



(26) 



where the dimensionless constants in the exponential are 
given by K = 12.0 and K\ — 5.1 for the system in 
question. 

We also tested how the properties of the system depend 
on the thickness of the wall separating the two media. For 
a fixed value of the final dielectric constant e-i — 2. we 
allowed d to vary. A sample of the results is shown in 
Fig. 10. As expected, for larger values of d the transition 
between the initial and final media is smoother, and the 
production rate is lower. As before, for d < 8 the final 
particle number depends exponentially on the expansion 
velocity, though for higher values of the wall thickness 
we observe deviations from this simple law. In the "thin 
wall" regime, the exponent 7 varies very little with d, as 
can be seen in the right plot in Fig. 10. For 1 < d < 7 
the change in 7 is less than 10%. 

Finally we checked whether these results depend on 
the specific shape of the wall profile. We replaced the sin 
wall by one where a linear function interpolates between 
the two dielectric constants. We observed that the results 
were virtually unchanged. 



V. CONCLUSIONS 

We have developed and tested the numerical tech- 
niques necessary to tackle the problem of particle pro- 
duction from a vacuum state in general scenarios. In 
particular, we have shown that the approach performs 
quite well when applied to cosmology, to the dynami- 
cal Casimir effect and to systems with complex optical 
properties. In some of these cases we were able to use the 
numerical data to obtain simple phcnomcnological laws. 
This illustrates the power of the method and shows how 
it can be used as a tool to approach problems for which 
there is little or no analytical information. The oppor- 
tunities are many, as was mentioned in the bulk of the 
paper. In the context of the dynamical Casimir effect, 
the first obvious step will be to generalise this approach 
to higher dimensions, where little is known analytically. 
This poses no obvious numerical difficulty, in contrast to 
the analytical problems raised by the loss of conformal 
invariance. Extensions to general types of geometry and 
the inclusion of a mass term in the theory should also 
be within reach of the method. Optical systems in 2D 
and 31? with alternative Lagrangians, and realistic phase 
transitions can also be subject to the same techniques. 
Finally, we should stress that since in general we are able 
to compute the full set of Bogolugov coefficients, we do 
not have to restrict the initial state to be in the vacuum. 
Time evolution of multi-particle states such as thermal 



configurations, can be easily simulated within this frame- 
work. 
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APPENDIX A: NUMERICAL ALGORITHM FOR 
MOVING BOUNDARY CONDITIONS 

In all simulations discussed above the scalar field wave 
equation was solved using a straight-forward Leap Frog 
algorithm (see for example [23]). Only the case of the 
moving mirror (Section III) requires more attention. For 
a box with time-dependent length, if the lattice spacing 
is kept constant, the number of points in the simulation 
domain must change and this may lead to large numerical 
errors. 

The Leap Frog (LF) discretization of the wave equation 
Eq. (15) is 



iL,„+i/2 = n ijTl _i/ 2 + dt v 2 

<fii,n+l = <Pi,n + dtHi , n +l/2j 



(Al) 



where i is the spatial lattice index and n the time step. 
As usual with LF-type algorithms, the field and its mo- 
mentum are defined respectively on integer and half- 
integer time-steps. As before the coordinate of the mov- 
ing boundary will be denoted by X(t), and at each time- 
step its value will not necessarily lie exactly on one of the 
lattice spatial points. This discrepancy will have to be 
taken into account. In particular, as X(t) increases new 
points will be introduced in the lattice and one must find 
a rule for assigning a field or momentum value to these. 

Let us start with the case where the new lattice point 
i + 1 lies in a half- integer time slice n + 1/2, as depicted 
in Fig. 11. 
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FIG. 11. Update algorithm for the moving boundary prob- 
lem: integer i denotes spatial lattice index and n is the 
time-step. For higher accuracy, we define auxiliary sites in 
both space and time directions as the boundary's trajectory 
- diagonal line - crosses the lattice. The separation between 
these and their adjacent integer lattice coordinates is denoted 
respectively, by a and (3. 

We start by differentiating in respect to time the 
boundary condition (f>(X(t),t) = 0, 



fa \{x(t),t) X{t) 



U(X(t),t) = 



(A2) 



This allows us to calculate the momentum at the bound- 
ary in terms of the spatial derivative of the field. Since 
the field vanishes at the boundary, we can approximate 
d x (p at the point (i+1, n+1/2) by — (<fii,n + 4 l i,n+i)/ (2 dx), 
using the fact that </>i jn +i can be evaluated in advance. 
This leads to: 



n 



x 



i+l, ra+l/2 



n+1/2 

2dx 



(<Pi,n+l + <t>i,n) 



(A3) 



Note however that here we implicitly assume that the 
boundary crosses the lattice at precisely the point (i + 
l,n + 1/2), which is not necessarily true. The above 
estimate can be improved by introducing an auxiliary 
fictitious lattice point (i + 1, n+1/2 — /3), the exact loca- 
tion where the mirror's trajectory intersects the lattice. 
The time difference between the crossing and time-step 
n + 1/2 is given by (3dt and can be evaluated by inter- 
polating the boundary trajectory between steps n and 
n+1/2. At the crossing point we have 



n 



x. 



i+l,n+l/2-/3 — 



n+1/2 

2dx 



[(1-2/3) &, n+1 + (1 + 2/3) </>i 



(A4) 



Note that by setting (3 = we recover the previous esti- 
mate Eq. (A3). Using Eq. (A4) we can then update the 
field at time-step n + 1 as 



dt (1/2 + /3) n i+li „ +1/2 _ /J 



(A5) 



The next momentum up-date, at n + 3/2, will be calcu- 
lated taking into account that boundary was placed at 
(i + l,n+ 1/2-/3). 

This method is easily generalised to the situation where 
the newly introduced lattice point lies on a integer time 
slice n + 1. We proceed in a similar way, defining an 
auxiliary momentum point at (i + 1, n+1/2 — /3), though 
in this case (3 will be negative. The following updates 
proceed exactly in the same way as described in Eqs. (A4) 
and (A5). 

When the trajectory of the mirror is such that no new 
points need to be introduced, we must still take some care 
when updating the fields near the boundary. In particu- 
lar, the calculation of the second derivative term V 2 cj> |j „ 



in Eq. (Al) should reflect the fact that the boundary may 
lie between lattice sites. Referring again to Fig. 11 we de- 
fine a as the distance in lattice units from point (i, n) to 
the mirror. Since the field vanishes at the boundary we 
will have 



aa 



-i,n - (1 + a)0»,r 



adx 2 



(A6) 



method 1 



By setting a = 1, this expression reduces to the less 
accurate estimate for the second derivative, where the 
boundary's position is identified with lattice point i + l. 
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FIG. 12. Numerical errors for three different methods for 
the first component of the momentum field IIi. A single mode 
ui — mr/Xo with n — 16 was evolved and compared with the 
corresponding analytical solution. The cavity expanded with 
constant velocity v = 0.9 from initial size Xo = 50. 

In order to check the accuracy of the algorithm, we 
have used the exact analytical solution Eq. (20) for the 
uniformly expanding box and compared it to the corre- 
sponding numerical solution. In Fig. 12 we show for a par- 
ticular solution, the numerical errors for three different 
methods based on the approximations discussed above. 
The error is defined as the (normalised) integral of the 
square of the difference between the exact and numeri- 
cal solutions. Algorithm 1 is obtained by setting (3 = 
and a = 1 in Eqs. (A4) and (A6). This corresponds 
to a straight forward implementation of the Leap Frog 
method with no smoothing at the boundary in cither time 
or space directions. Method 2 includes the corrections for 
the momentum at the boundary but non-improved spa- 
tial derivatives (a = 1). Finally, in method 3 both a and 
(3 are evaluated at each time-step according to the mir- 
rors's position. The improvement in accuracy is clear: 
algorithm 3 decreases the errors by a factor of 10 over 
method 2 and by a factor of 25 over method 1. For the 
purposes of the simulations in this work this is enough. 
Direct observation of the numerical solutions shows that 
the moving boundary introduces high-frequency discon- 
tinuities in the fields (in particular in the momentum). 
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These are considerably reduced by method 3, and in the 
end they tend to average out when the Bogolugov coef- 
ficients are evaluated, leading to reasonably small errors 
in the final physical quantities. 
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